第六章:多層次分析
鄭中平、許清芳 2024 三月 03
#整體設定,含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")
資料
dta <- read.csv("../Data/PIRLS2016TW.csv", stringsAsFactors = TRUE)
[1] 4112 8
|
學生
|
老師
|
學校
|
性別
|
書本數
|
女性比率
|
社經地位
|
閱讀成績
|
|
1
|
101
|
1
|
男
|
18
|
0.5
|
差不多
|
613.1
|
|
2
|
101
|
1
|
男
|
5
|
0.5
|
差不多
|
610.7
|
|
3
|
101
|
1
|
男
|
63
|
0.5
|
差不多
|
534.7
|
|
4
|
101
|
1
|
男
|
251
|
0.5
|
差不多
|
613.6
|
|
5
|
101
|
1
|
男
|
5
|
0.5
|
差不多
|
654.3
|
|
6
|
101
|
1
|
男
|
18
|
0.5
|
差不多
|
680.9
|
描述統計
#重新分派「社經地位」的水準順序
dta <- dta |>
dplyr::mutate(社經地位 = fct_relevel(社經地位,
c("弱勢多", "差不多", "富有多")))
#程式報表6.2
dta[,-c(1:3)] |>
gtsummary::tbl_summary(by = '社經地位',
statistic = list(all_continuous() ~ "{mean} ({sd})")) |>
gtsummary::add_overall()
| Characteristic |
Overall, N = 4,112 |
弱勢多, N = 466 |
差不多, N = 2,498 |
富有多, N = 1,148 |
| 性別 |
|
|
|
|
| 女 |
1,972 (48%) |
233 (50%) |
1,208 (48%) |
531 (46%) |
| 男 |
2,140 (52%) |
233 (50%) |
1,290 (52%) |
617 (54%) |
| 書本數 |
|
|
|
|
| 5 |
624 (15%) |
133 (29%) |
382 (15%) |
109 (9.5%) |
| 18 |
871 (21%) |
124 (27%) |
556 (22%) |
191 (17%) |
| 63 |
1,278 (31%) |
127 (27%) |
801 (32%) |
350 (30%) |
| 151 |
692 (17%) |
47 (10%) |
405 (16%) |
240 (21%) |
| 251 |
647 (16%) |
35 (7.5%) |
354 (14%) |
258 (22%) |
| 女性比率 |
0.48 (0.07) |
0.50 (0.08) |
0.48 (0.06) |
0.46 (0.07) |
| 閱讀成績 |
564 (62) |
542 (62) |
561 (62) |
579 (59) |
生態相關(Ecological
correlation)
#看看書本數, 閱讀成績相關,以學生層次分數計算
dta |> dplyr::select(書本數, 閱讀成績) |> cor()
|
|
書本數
|
閱讀成績
|
|
書本數
|
1.0000
|
0.2731
|
|
閱讀成績
|
0.2731
|
1.0000
|
#計算各校書本數與閱讀成績平均分數
dta_m <- dta |>
dplyr::group_by(學校) |>
dplyr::summarize(書本數平均 = mean(書本數),
閱讀成績平均 = mean(閱讀成績))
#以學校層次分數計算生態相關(ecological correlation)
dta_m |> dplyr::select(書本數平均, 閱讀成績平均) |> cor()
|
|
書本數平均
|
閱讀成績平均
|
|
書本數平均
|
1.0000
|
0.6444
|
|
閱讀成績平均
|
0.6444
|
1.0000
|
繪圖
#直方圖可以瞭解分數的變異
#這邊以學生成績畫直方圖,看到學生間的變異
#再以學生平均(各校)畫直方圖,看到校間變異
#圖6.1
ggplot() +
stat_bin(data=dta,
aes(閱讀成績, after_stat(density)),
alpha=.2,
color='white',
binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)))+
stat_bin(data=dta_m,
aes(閱讀成績平均, after_stat(density)),
alpha=.2,
color='black',
binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)))+
labs(x = '閱讀成績',
y = '密度',
title = '學生閱讀成績, 各校閱讀成績平均分數')
#繪製學校內學生分數標準差的直方圖,可以知道不同校內變異大小的學校各佔多少
#圖6.2
dta |> dplyr::group_by(學校) |>
dplyr::summarize(閱讀成績標準差 = sd(閱讀成績)) |>
ggplot()+
stat_bin(aes(閱讀成績標準差, after_stat(density)),
color='black',
fill='gray',
binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)))+
labs(x = '學校閱讀成績標準差',
y = '密度')
#此圖縱軸是學校,橫軸是學校成績的平均與標準誤,並且依平均由大到小排列學校
#圖6.3
dta |>
ggplot()+
aes(y = 閱讀成績, x = reorder(學校, 閱讀成績, mean)) +
stat_summary(fun.data = "mean_se", size=rel(.1), alpha=.2) +
geom_hline(yintercept=mean(dta$閱讀成績),
linetype='dotted',
col=8)+
coord_flip() +
labs(y = '閱讀成績平均分數',
x = '學校')+
theme(axis.text.y=element_blank())
#預計要先看學校社經地位對閱讀成績的影響
#但學校社經地位為「富有多」時,女性比率較低,因此先畫圖呈現此一現象
#圖6.4
ggplot(dta,
aes(x=女性比率, after_stat(density), color=社經地位)) +
geom_density()+
scale_color_grey()+
scale_x_continuous(breaks=seq(.2,.8,by=.1))+
labs(x="教師班女性比率",
y="密度")+
theme(legend.position='top')
#依社經地位,各教師班女性比率與閱讀成績的散佈圖,並加上迴歸線
#圖6.5
ggplot(dta,
aes(女性比率, 閱讀成績)) +
geom_point(size=rel(.5), alpha=.5)+
stat_smooth(method='lm',
formula = y ~ x,
se =FALSE,
linewidth=.5,
col=8)+
scale_color_grey()+
facet_wrap(vars(社經地位))+
labs(x="教師班女性比率",
y="閱讀成績")
#書本數與閱讀成績的散佈圖,加上以校為單位的迴歸線
#圖6.6
ggplot(dta,
aes(書本數, 閱讀成績, group=學校)) +
geom_point(size=rel(.2), alpha=.2)+
stat_smooth(method='lm',
formula = y ~ x,
se =FALSE,
linewidth=.2,
col=8,
alpha=.1)+
scale_x_continuous(trans='log2')+
labs(x="書籍數目",
y="閱讀成績")
#看一下各班性別比率的統計
with(dta, 女性比率) |> summary()
|
Min.
|
1st Qu.
|
Median
|
Mean
|
3rd Qu.
|
Max.
|
|
0.25
|
0.4286
|
0.4783
|
0.4796
|
0.5185
|
0.7917
|
#各校以書本數預測閱讀成績數時的截距與斜率,分各校畫圖
#女性比率可能是兩者關係混淆變項,此處選擇女性比率在中位數附近者
#圖6.7
dta |> dplyr::filter(女性比率 > .455, 女性比率 < .48) |>
ggplot() +
aes(x = 書本數, y = 閱讀成績) +
geom_point(size = rel(.6)) +
stat_smooth(method = 'lm', formula = y ~ x, se = F,
linewidth=.5, col=1) +
facet_wrap(vars(學校))+
scale_x_continuous(trans = 'log2')+
labs(x="書籍數目",
y="閱讀成績",
subtitle='學校')+
theme(legend.position='top')
#區分社經地位與性別,以蜂群圖(beeswarm)繪製閱讀成績
#圖6.8
ggplot(data=dta,
aes(x=性別, y=閱讀成績))+
geom_boxplot()+
ggbeeswarm::geom_beeswarm(size=rel(.2)) +
facet_wrap(vars(社經地位))
#區分社經地位與性別,繪製書本數與閱讀成績關係
#圖6.9
ggplot(data=dta,
aes(x=書本數, y=閱讀成績))+
stat_smooth(method='lm',
formula = y ~ x,
col='black')+
stat_summary(fun.data="mean_se", size=rel(.5))+
facet_grid(性別 ~ 社經地位)+
scale_x_continuous(trans='log2')
多層次分析
#分析時請確認,教師與學校編號,都已經被視為 factor
#在此將書本數置中
dta <- dta |>
dplyr::mutate(老師 = as.factor(老師),
學校 = as.factor(學校),
書本數置中 = log2(書本數) |> scale(scale=FALSE) |>
as.vector())
#只包含隨機截距的模型
#程式報表6.3
m0 <- lme4::lmer(閱讀成績 ~ 1 + ( 1 | 學校/老師 ),
data = dta,
na.action=na.omit)
summary(m0, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: 閱讀成績 ~ 1 + (1 | 學校/老師)
Data: dta
REML criterion at convergence: 45398
Scaled residuals:
Min 1Q Median 3Q Max
-4.19 -0.59 0.08 0.66 3.10
Random effects:
Groups Name Variance Std.Dev.
老師:學校 (Intercept) 97 9.85
學校 (Intercept) 259 16.10
Residual 3479 58.98
Number of obs: 4112, groups: 老師:學校, 171; 學校, 145
Fixed effects:
Estimate Std. Error t value
(Intercept) 561.63 1.83 308
#同一個模型的另一種寫法
#雖未直接設定老師與學校的巢套關係,但資料中仍可以看出
#建議使用這種寫法
m0 <- lme4::lmer(閱讀成績 ~ 1 + ( 1 | 老師 ) + ( 1 | 學校 ),
data = dta,
na.action=na.omit)
summary(m0)
#參數估計
#程式報表6.4
m0 |> broom::tidy()
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
561.626
|
1.826
|
307.5
|
|
ran_pars
|
老師:學校
|
sd\_\_(Intercept)
|
9.847
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
16.099
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
58.981
|
NA
|
NA
|
加入學校層級解釋變項
#加入社經地位(學校層級獨變項)
#程式報表6.5
broom::tidy(m1 <- update(m0, . ~ . + 社經地位))
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
539.994
|
4.276
|
126.297
|
|
fixed
|
NA
|
社經地位差不多
|
20.650
|
4.729
|
4.366
|
|
fixed
|
NA
|
社經地位富有多
|
37.082
|
5.304
|
6.992
|
|
ran_pars
|
老師:學校
|
sd\_\_(Intercept)
|
9.588
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
12.072
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
58.976
|
NA
|
NA
|
加入教師班層級解釋變項
##加入女性比率(教師班層級)及其與社經地位交互作用(跨層次交互作用)
#程式報表6.6
broom::tidy(m2 <- update(m1, . ~ . + 女性比率 + 女性比率:社經地位))
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
556.981
|
24.82
|
22.4407
|
|
fixed
|
NA
|
社經地位差不多
|
-11.780
|
28.70
|
-0.4104
|
|
fixed
|
NA
|
社經地位富有多
|
12.676
|
30.93
|
0.4098
|
|
fixed
|
NA
|
女性比率
|
-33.638
|
48.32
|
-0.6961
|
|
fixed
|
NA
|
社經地位差不多:女性比率
|
65.472
|
56.59
|
1.1569
|
|
fixed
|
NA
|
社經地位富有多:女性比率
|
49.472
|
62.06
|
0.7971
|
|
ran_pars
|
老師:學校
|
sd\_\_(Intercept)
|
9.496
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
12.378
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
58.969
|
NA
|
NA
|
加入學生層級解釋變項
##加入書本數變項(學生層級、置中)、性別(學生層級)及其與社經地位交互作用(跨層次交互作用)
#程式報表6.7
broom::tidy(m3 <- update(m2, . ~ . + 性別 + 書本數置中 + 性別:社經地位))
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
575.054
|
22.1454
|
25.9672
|
|
fixed
|
NA
|
社經地位差不多
|
-22.801
|
25.5041
|
-0.8940
|
|
fixed
|
NA
|
社經地位富有多
|
-4.054
|
27.3719
|
-0.1481
|
|
fixed
|
NA
|
女性比率
|
-41.719
|
42.2983
|
-0.9863
|
|
fixed
|
NA
|
性別男
|
-10.791
|
5.3239
|
-2.0270
|
|
fixed
|
NA
|
書本數置中
|
9.679
|
0.5006
|
19.3337
|
|
fixed
|
NA
|
社經地位差不多:女性比率
|
69.482
|
49.4078
|
1.4063
|
|
fixed
|
NA
|
社經地位富有多:女性比率
|
51.996
|
53.8922
|
0.9648
|
|
fixed
|
NA
|
社經地位差不多:性別男
|
3.079
|
5.7938
|
0.5315
|
|
fixed
|
NA
|
社經地位富有多:性別男
|
5.802
|
6.3121
|
0.9192
|
|
ran_pars
|
老師:學校
|
sd\_\_(Intercept)
|
7.786
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
9.326
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
56.649
|
NA
|
NA
|
模型比較
#比較四個模型
#程式報表6.8
anova(m0, m1, m2, m3)
refitting model(s) with ML (instead of REML)
|
|
npar
|
AIC
|
BIC
|
logLik
|
deviance
|
Chisq
|
Df
|
Pr(\>Chisq)
|
|
m0
|
4
|
45409
|
45434
|
-22700
|
45401
|
NA
|
NA
|
NA
|
|
m1
|
6
|
45369
|
45407
|
-22678
|
45357
|
43.948
|
2
|
0.0000
|
|
m2
|
9
|
45373
|
45430
|
-22677
|
45355
|
1.823
|
3
|
0.6099
|
|
m3
|
13
|
45008
|
45090
|
-22491
|
44982
|
373.097
|
4
|
0.0000
|
#看一下m3模型表現(AIC, BIC, conditional R2, marginal R2)
#程式報表6.9前
performance::performance(m3) |> dplyr::select(-c(2,6:8))
|
AIC
|
BIC
|
R2_conditional
|
R2_marginal
|
|
44957
|
45039
|
0.156
|
0.1172
|
組內相關(Intraclass
correlations)
#比較m0與m3的ICC,看看學生間的相似性有多少被解釋變項解釋掉了
#程式報表6.9後
performance::icc(m0)
|
ICC_adjusted
|
ICC_unadjusted
|
optional
|
|
0.0929
|
0.0929
|
FALSE
|
|
ICC_adjusted
|
ICC_unadjusted
|
optional
|
|
0.044
|
0.0388
|
FALSE
|
#抽取變異成分,計算教師與學校可以解釋部分(ICCs)
print(vc <- lme4::VarCorr(m3) |> as.data.frame(), comp = 'Variance')
grp var1 var2 vcov sdcor
1 老師:學校 (Intercept) <NA> 60.63 7.786
2 學校 (Intercept) <NA> 86.97 9.326
3 Residual <NA> <NA> 3209.09 56.649
#換成百分比
vc["vcov"]/sum(vc["vcov"])
|
vcov
|
|
0.0181
|
|
0.0259
|
|
0.9560
|
以圖形顯示模型效果
#圖示固定效果
#將截距去除,圖較更易懂
#圖6.10左
p1 <-
broom::tidy(m3, conf.int = TRUE) |>
dplyr::filter(effect=='fixed') |>
dplyr::slice(-1) |>
ggplot() +
aes(x=estimate, y=term, xmin = conf.low, xmax = conf.high, height = 0) +
geom_point(size=rel(3)) +
geom_vline(xintercept = 0, linetype="dashed") +
geom_errorbarh() +
labs(x = '估計值和信賴區間',
y = '模型固定效果變項(去除截距)',
subtitle = '依變項:閱讀成績')
#也可以用以下程式碼得出類似圖形
#sjPlot::plot_model(m3)
#圖示隨機效果(估計值與信賴區間)
#圖6.10右
p2<-
lme4::ranef(m3) |> as.data.frame() |>
ggplot() +
aes(y=grp, x=condval, color=grpvar) +
geom_point(size=rel(.8)) +
facet_wrap(vars(term), scale='free') +
geom_errorbarh(aes(xmin=condval-2*condsd,
xmax=condval+2*condsd),
linewidth=rel(.4), alpha=.2)+
scale_color_grey()+
labs(x="隨機效果的條件平均數",
y="")+
guides(color=guide_legend(title="單位"))+
theme(axis.text.y=element_blank(),
legend.position="top")
p1+p2
模型診斷(檢驗假設)
殘差
#圖6.11
gglm::gglm(m3, theme = ggplot2::theme_minimal())
隨機效果QQ圖
#圖6.12
lme4::ranef(m3) |> as.data.frame() |>
ggplot() +
aes(sample=condval)+
stat_qq(size=rel(.5)) +
stat_qq_line() +
facet_wrap(vars(grpvar)) +
labs(x="常態分位數",
y="隨機效果的條件平均數")
模型配適
#比對社經地位與女性比率與閱讀成績關係的實際資料與預測值
#圖6.13
ef_x <- effects::effect(term=c("女性比率", "社經地位"), mod=m3) |>
as.data.frame()
ggplot() +
geom_point(data=dta, aes(女性比率, 閱讀成績), alpha=.2, pch=1) +
geom_point(data=ef_x, aes(x=女性比率, y=fit)) +
geom_ribbon(data=ef_x, aes(x=女性比率,
ymin=lower,
ymax=upper), fill=8) +
geom_line(data=ef_x, aes(x=女性比率, y=fit)) +
facet_wrap(vars(社經地位)) +
labs(x="女性比率",
y="閱讀成績")
#比對書本數與女性比率與閱讀成績關係的實際資料與預測值
#圖6.14
broom::augment(m3, dta) |>
dplyr::filter(女性比率 > .455, 女性比率 < .48) |>
ggplot() +
aes(x = 書本數, y = 閱讀成績) +
geom_point(size = rel(.5)) +
stat_smooth(method = 'lm',
formula = y ~ x,
se = F, col=1,
linewidth=.3) +
stat_smooth(aes(y=.fitted),
method = 'lm',
formula = y ~ x,
se = F, col=2,
linewidth=.3) +
facet_wrap(vars(學校)) +
scale_x_continuous(trans='log2') +
labs(x="書本數",
y="閱讀成績",
subtitle='學校')+
theme(legend.position='top')
複雜的多層次模型
#隨機斜率模型,不考慮截距與斜率相關
#程式報表6.10
broom::tidy(m4 <- update(m3, . ~ . + (書本數置中 - 1 | 學校)))
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
575.117
|
22.1793
|
25.9304
|
|
fixed
|
NA
|
社經地位差不多
|
-22.175
|
25.5396
|
-0.8683
|
|
fixed
|
NA
|
社經地位富有多
|
-4.495
|
27.3870
|
-0.1641
|
|
fixed
|
NA
|
女性比率
|
-42.124
|
42.3770
|
-0.9940
|
|
fixed
|
NA
|
性別男
|
-10.873
|
5.3168
|
-2.0451
|
|
fixed
|
NA
|
書本數置中
|
9.764
|
0.5606
|
17.4177
|
|
fixed
|
NA
|
社經地位差不多:女性比率
|
68.427
|
49.4958
|
1.3825
|
|
fixed
|
NA
|
社經地位富有多:女性比率
|
53.151
|
53.9190
|
0.9858
|
|
fixed
|
NA
|
社經地位差不多:性別男
|
3.285
|
5.7864
|
0.5678
|
|
fixed
|
NA
|
社經地位富有多:性別男
|
5.930
|
6.3023
|
0.9410
|
|
ran_pars
|
老師.學校
|
sd\_\_(Intercept)
|
7.836
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
9.133
|
NA
|
NA
|
|
ran_pars
|
學校.1
|
sd\_\_書本數置中
|
2.886
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
56.422
|
NA
|
NA
|
#原始報表
summary(m4, corr=FALSE)
Linear mixed model fit by REML ['lmerMod']
Formula: 閱讀成績 ~ (1 | 學校/老師) + 社經地位 + 女性比率 +
性別 + 書本數置中 + (書本數置中 - 1 | 學校) +
社經地位:女性比率 + 社經地位:性別
Data: dta
REML criterion at convergence: 44927
Scaled residuals:
Min 1Q Median 3Q Max
-4.682 -0.595 0.067 0.668 3.351
Random effects:
Groups Name Variance Std.Dev.
老師.學校 (Intercept) 61.41 7.84
學校 (Intercept) 83.41 9.13
學校.1 書本數置中 8.33 2.89
Residual 3183.43 56.42
Number of obs: 4112, groups: 老師:學校, 171; 學校, 145
Fixed effects:
Estimate Std. Error t value
(Intercept) 575.117 22.179 25.93
社經地位差不多 -22.175 25.540 -0.87
社經地位富有多 -4.495 27.387 -0.16
女性比率 -42.124 42.377 -0.99
性別男 -10.873 5.317 -2.05
書本數置中 9.764 0.561 17.42
社經地位差不多:女性比率 68.427 49.496 1.38
社經地位富有多:女性比率 53.151 53.919 0.99
社經地位差不多:性別男 3.285 5.786 0.57
社經地位富有多:性別男 5.930 6.302 0.94
#隨機斜率模型,考慮截距與斜率相關
#程式報表6.11
broom::tidy(m5 <- lme4::lmer(閱讀成績 ~ 書本數置中 +
(1 | 老師) +
(1 + 書本數置中 | 學校),
data = dta,
na.action = na.omit))
|
effect
|
group
|
term
|
estimate
|
std.error
|
statistic
|
|
fixed
|
NA
|
(Intercept)
|
562.7550
|
1.4559
|
386.54
|
|
fixed
|
NA
|
書本數置中
|
10.0581
|
0.5605
|
17.94
|
|
ran_pars
|
老師
|
sd\_\_(Intercept)
|
7.8414
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_(Intercept)
|
11.1882
|
NA
|
NA
|
|
ran_pars
|
學校
|
cor\_\_(Intercept).書本數置中
|
0.1232
|
NA
|
NA
|
|
ran_pars
|
學校
|
sd\_\_書本數置中
|
2.9294
|
NA
|
NA
|
|
ran_pars
|
Residual
|
sd\_\_Observation
|
56.5503
|
NA
|
NA
|
#呈現此時的 ICC
performance::icc(m5)
|
ICC_adjusted
|
ICC_unadjusted
|
optional
|
|
0.0633
|
0.0575
|
FALSE
|
縮減(Shrinkage)
#比較資料、迴歸、多層次迴歸間截距與斜率的關係
#收集資料準備繪圖
betas <- lme4::lmList(閱讀成績 ~ 書本數置中 | 學校, data=dta) |>
coef() |>
as_tibble( rownames = "rowname")
bs <- m5 |> coef() |>
purrr::pluck("學校") |>
as_tibble( rownames = "rowname")
b_df <- dplyr::inner_join(betas, bs, by='rowname')
names(b_df) <- c("id", "beta0", "beta1", "b0", "b1")
#圖6.15
ggplot(data=b_df,
aes(beta0, beta1)) +
geom_point(aes(x=lme4::fixef(m5)[1],
y=lme4::fixef(m5)[2]),
col=1, size=rel(2))+
geom_point(size=rel(1.5), shape=21, col=8) +
geom_point(aes(b0, b1),
size=rel(1), shape=16, col=1) +
geom_segment(aes(x=beta0, y=beta1,
xend=b0, yend=b1),
alpha=.2, size=rel(.5),
arrow=arrow(angle=15,
length=unit(.2, "cm")))+
stat_ellipse(col=8, linetype=3)+
stat_ellipse(aes(b0, b1), col=1, linetype=1)+
labs(x="Intercept estimats",
y="Slope estimates")
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
##附錄:多層次分析與隨機效果ANOVA 間的關聯
#程式報表6.12
aov(閱讀成績 ~ Error(老師), data=dta) |> summary()
Error: 老師
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 170 2063990 12141
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3941 13700482 3476
aggregate(學生 ~ 老師, data=dta, FUN=length) |> pluck('學生') |> mean()
[1] 24.05
[1] 360.3
References
紀馥安, 許清芳(2012).運用開放軟體 R
處理大型教育資料庫.當代教育研究季刊, 23(4), 121-153.
結束
顯示演練單元信息
R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)
Matrix products: default
locale:
[1] LC_COLLATE=Chinese (Traditional)_Taiwan.utf8
[2] LC_CTYPE=Chinese (Traditional)_Taiwan.utf8
[3] LC_MONETARY=Chinese (Traditional)_Taiwan.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=Chinese (Traditional)_Taiwan.utf8
time zone: Asia/Taipei
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] brolgar_1.0.0.9000 ggmirt_0.1.0 mlogit_1.1-1
[4] dfidx_0.0-5 dispmod_1.2 metaviz_0.3.1
[7] metafor_4.4-0 numDeriv_2016.8-1.1 metadat_1.2-0
[10] ngramr_1.9.3 rentrez_1.2.3 psych_2.3.9
[13] psychometric_2.4 multilevel_2.7 nlme_3.1-163
[16] cowplot_1.1.1 eRm_1.0-4 qgraph_1.9.8
[19] semTools_0.5-6 CTT_2.3.3 tidySEM_0.2.6
[22] OpenMx_2.21.11 semPlot_1.1.6 lavaanPlot_0.6.2
[25] lavaan_0.6-16 effects_4.2-2 carData_3.0-5
[28] lme4_1.1-35.1 broom.mixed_0.2.9.4 plyr_1.8.9
[31] mediation_4.5.0 sandwich_3.1-0 mvtnorm_1.2-4
[34] Matrix_1.6-4 MBESS_4.9.3 vctrs_0.6.2
[37] rsvg_2.6.0 DiagrammeRsvg_0.1 DiagrammeR_1.0.10
[40] bda_17.1.2 boot_1.3-28.1 insight_0.19.7
[43] lm.beta_1.7-2 see_0.8.1 performance_0.10.8
[46] ggokabeito_0.1.0 ggrepel_0.9.4 gglm_1.0.2
[49] ggpubr_0.6.0 gghighlight_0.4.0 ggridges_0.5.4
[52] ggbeeswarm_0.7.2 KernSmooth_2.23-22 ggExtra_0.10.1
[55] ggeffects_1.3.2 webshot2_0.1.1 patchwork_1.1.3
[58] ragg_1.2.7 GGally_2.2.0 moments_0.14.1
[61] corrplot_0.92 cocor_1.1-4 parameters_0.21.3
[64] Hmisc_5.1-1 confintr_1.0.2 kableExtra_1.3.4
[67] MASS_7.3-60 statmod_1.5.0 flextable_0.9.4
[70] printr_0.3 here_1.0.1 sjPlot_2.8.15
[73] broom_1.0.5 gtsummary_1.7.2 jtools_2.2.2
[76] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
[79] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
[82] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4
[85] tidyverse_2.0.0 remotes_2.4.2.1 devtools_2.4.5
[88] usethis_2.2.2 pacman_0.5.1
loaded via a namespace (and not attached):
[1] sjmisc_2.8.9 urlchecker_1.0.1 kutils_1.73
[4] nnet_7.3-19 rstan_2.32.3 bain_0.2.10
[7] digest_0.6.31 png_0.1-8 corpcor_1.6.10
[10] bayestestR_0.13.1 httpcode_0.3.0 parallelly_1.36.0
[13] permute_0.9-7 fontLiberation_0.1.0 reshape2_1.4.4
[16] httpuv_1.6.13 withr_2.5.2 xfun_0.39
[19] ellipsis_0.3.2 survival_3.5-7 commonmark_1.9.0
[22] memoise_2.0.1 crul_1.4.0 emmeans_1.8.9
[25] profvis_0.3.8 systemfonts_1.0.5 zoo_1.8-12
[28] MplusAutomation_1.1.0 gtools_3.9.5 V8_4.4.1
[31] pbapply_1.7-2 Formula_1.2-5 datawizard_0.9.0
[34] promises_1.2.1 httr_1.4.7 rstatix_0.7.2
[37] rockchalk_1.8.157 globals_0.16.2 ps_1.7.5
[40] rstudioapi_0.15.0 miniUI_0.1.1.1 generics_0.1.3
[43] base64enc_0.1-3 processx_3.8.3 curl_5.2.0
[46] mitools_2.4 quadprog_1.5-8 xtable_1.8-4
[49] evaluate_0.23 hms_1.1.3 colorspace_2.1-0
[52] visNetwork_2.1.2 lmtest_0.9-40 magrittr_2.0.3
[55] later_1.3.2 lattice_0.21-9 tsibble_1.1.3
[58] future.apply_1.11.0 matrixStats_1.2.0 XML_3.99-0.16
[61] survey_4.2-1 texreg_1.39.3 pillar_1.9.0
[64] StanHeaders_2.26.28 compiler_4.3.2 stringi_1.7.12
[67] minqa_1.2.6 crayon_1.5.2 abind_1.4-5
[70] mathjaxr_1.6-0 modelr_0.1.11 chromote_0.1.2
[73] codetools_0.2-19 textshaping_0.3.7 nonnest2_0.5-6
[76] openssl_2.1.1 QuickJSR_1.0.8 mime_0.12
[79] anytime_0.3.9 splines_4.3.2 markdown_1.12
[82] Rcpp_1.0.11 fastDummies_1.7.3 knitr_1.45
[85] utf8_1.2.3 pbivnorm_0.6.0 fs_1.6.3
[88] listenv_0.9.0 checkmate_2.2.0 Rdpack_2.6
[91] pkgbuild_1.4.3 openxlsx_4.2.5.2 ggsignif_0.6.4
[94] estimability_1.4.1 tzdb_0.4.0 lpSolve_5.6.20
[97] svglite_2.1.3 bayesplot_1.10.0 pkgconfig_2.0.3
[100] tools_4.3.2 cachem_1.0.8 rbibutils_2.2.16
[103] viridisLite_0.4.2 blavaan_0.5-2 rvest_1.0.3
[106] DBI_1.1.3 fastmap_1.1.1 rmarkdown_2.25
[109] scales_1.3.0 grid_4.3.2 gt_0.10.0
[112] sass_0.4.8 officer_0.6.3 coda_0.19-4.1
[115] ggstats_0.5.1 tmvnsim_1.0-2 RANN_2.6.1
[118] rpart_4.1.21 farver_2.1.1 mgcv_1.9-0
[121] gsubfn_0.7 yaml_2.3.7 foreign_0.8-85
[124] sjstats_0.18.2 cli_3.6.1 stats4_4.3.2
[127] webshot_0.5.5 dbscan_1.1-12 lifecycle_1.0.4
[130] fabletools_0.3.4 askpass_1.2.0 sessioninfo_1.2.2
[133] backports_1.4.1 timechange_0.2.0 gtable_0.3.4
[136] progressr_0.14.0 parallel_4.3.2 jsonlite_1.8.8
[139] vegan_2.6-4 glasso_1.11 proto_1.0.0
[142] zip_2.3.0 RcppParallel_5.1.7 GPArotation_2023.11-1
[145] highr_0.10 sem_3.1-15 loo_2.6.0
[148] distributional_0.3.2 dcurver_0.9.2 pander_0.6.5
[151] shiny_1.8.0 lisrelToR_0.1.5 htmltools_0.5.7
[154] sjlabelled_1.2.0 glue_1.6.2 gfonts_0.2.0
[157] broom.helpers_1.14.0 gdtools_0.3.5 rprojroot_2.0.4
[160] mirt_1.41 mnormt_2.1.1 jpeg_0.1-10
[163] gridExtra_2.3 igraph_1.6.0 R6_2.5.1
[166] arm_1.13-1 fdrtool_1.2.17 labeling_0.4.3
[169] Deriv_4.1.3 CompQuadForm_1.4.3 cluster_2.1.4
[172] pkgload_1.3.3 nloptr_2.0.3 rstantools_2.3.1.1
[175] tidyselect_1.2.0 vipor_0.4.5 htmlTable_2.4.2
[178] xml2_1.3.4 inline_0.3.19 fontBitstreamVera_0.1.1
[181] car_3.1-2 future_1.33.0 munsell_0.5.0
[184] mi_1.1 furrr_0.3.1 fontquiver_0.2.1
[187] data.table_1.14.8 websocket_1.4.1 htmlwidgets_1.6.4
[190] RColorBrewer_1.1-3 rlang_1.1.1 uuid_1.1-1
[193] fansi_1.0.4 beeswarm_0.4.0